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ELF/VLF Electromagnetic Detection and Characterization 
of Deeply Buried Targets 


This is the Final Report for the SBIR Phase II effort entitled " ELF/VLF Electromagnetic 
Detection and Characterization of Deeply Buried Targets" sponsored by the Defense 
Advanced Research Projects Agency. The various electromagnetic calculational tools 
coded as MATLAB routines delivered with the report constitute the goal of our effort. 
The report is in five parts: (1) a description of the Sommerfeld routines; (2) a description 
of the menu interface routines; (3) a description of the Hill / King buried wire routines; 
(4) an appendix of the Sommerfeld integrals used; (5) an appendix describing the CD- 
ROM containing the MATLAB codes. 

Sommerfeld Calculation Routines for a Two-Laver Earth 

The configuration for the Sommerfeld routines is sketched below: 


z 



To calculate the electric or magnetic fields due to a dipole at a set of points, 
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MATLAB function routines are used: getElecDipoleField or getMagDipoleField; 
these functions are called in MATLAB as: 


[EField,HField] = ... 

getElecDipoleField(EFieldOn,HFieldOn,ILdipMom,. . . 
RSourcePoint,REvalPoints,f rec lVal, . . . 

zLayer1,eDiel eel,Conduct1,zLayer2,eDielec2,C° n duct2); 


[EField,HField] = ... 

getMagDipoleField(EFieldOn,HFieldOn,IAdipMom,.. . 
RSourcePoint,REvalPoints,freqVal,... 

zLayer1,eDiel eel,Conductl,zLayer2,eDielec2,C° n duct2); 

The description for most of the arguments is the same for both routines; the routines' 
input arguments differ only in the definition of the dipole moments. 


Below are the input 
EFieldOn 

HFieldOn 

RSourcePoint 


REvalPoints 


freqVal 


argument descriptions : 

- Logical variable: ON (=1) -> Calculate the E Field; 

OFF (=0) -> Don't Calculate the E Field 

- Logical variable: ON (=1) -> Calculate the H Field; 

OFF (=0) -> Don't Calculate the H Field 

- Location of the dipole source in meters for user defined 
Cartesian coordinates; RSourcePoint(l / 3), 

i.e. RSourcePoint = [xsrc,ysrc,zsrc] 

- Points to evaluate the E and / or H Fields in meters; uses the same 
(user defined) Cartesian coordinate system as RSourcePoint 
(above); REvalPoints(l:N,3), where N is the number of eval points 
i.e. REvalPoints(l,:) = [xEv(l) / yEv(l),zEv(l)] / 

REvalPoints(2 / :) = [xEv(2) / yEv(2),zEv(2)] / 

REvalPoints(3,:) = [xEv(3),yEv(3),zEv(3)], 

.etc. 

N.B. For a single call of either routine the REvalPoints selected 
must all be in a single medium; if the E/H Field values are desired 
in both Air and Region l, this may be accomplished in two 
sequential calls. Also note: Fields cannot be calculated in Region 2, 
as this is intended to model a substrate below the air and ground 
regions of sources and observation points 

- Frequency in Hz 
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zLayerl 

eDielecl 

Conductl 

zLayer2 

eDielec2 

Conduct2 

ILdipMom 

IAdipMom 


- z value [meters] of Air <--> Region 1 interface in the user defined 
Cartesian coordinate system 

- Relative dielectric constant of Region 1 medium 

- Conductivity of Region 1 medium [mhos / m] 

- z value [meters] of Region 1 <—> Region 2 interface in the user 
defined Cartesian coordinate system 

- Relative dielectric constant of Region 2 medium 

- Conductivity of Region 2 medium [mhos / m] 

- Electric dipole moment [amp-m]; used only for the 
function getElecDipoleField 

- Magnetic dipole moment [amp-m 2 ]; used only for the 
function getMagDipoleField 


Below are the 
EField 


output value descriptions : 

- The value of the E Field due to the dipole at the specified 
evaluation points (REvalPoints) in volt/m. The E Field is given in 
Cartesian coordinates. So, 

EField(l / :) = [Ex,Ey,Ez] at REvalPointsCl,:) 

EField(2 / :) = [Ex,Ey,Ez] at REvalPoints(2 / :) 

EField(3 / :) = [Ex,Ey,Ez] at REvalPoints(3 / :) 

.etc. 


FIField - The value of the H Field due to the dipole at the specified 

evaluation points (REvalPoints) in amp/m. The H Field is given in 
Cartesian coordinates. The array structure and correspondence 
between FIField and REvalPoints is the same as for EField above. 


There are a number of MATFAB functions required by getElecDipoleField and 
getMagDipoleField; these functions are supplied on the CD-ROM and they must be 
kept in the same folder as getElecDipoleField and getMagDipoleField. 

Example: Simple Buried HED Calculation 

As a first example we give below the code for calculation of an HED in Region 1 at a 
depth of 0.5 meter below the surface at 2 kHz. The relative dielectric constant of Region 
1 is 10 and the conductivity is 0.01 mhos/m; the relative dielectric constant of Region 2 
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is 20 and the conductivity is 0.05 mhos/m. The interface between the Air and Region 1 
corresponds to z = 0 m; the interface between the Region 1 and Region 2 corresponds 
to z = -100 m. We assume the electric dipole strength is 0.5 amp-m. The dipole is 
centered at x = 0 m, y = 0 m, and z = -0.5 m in our coordinate system and the dipole 
vector is along the y-axis. The electric and magnetic fields are evaluated on a square 
grid 30 m x 30 m a distance of 1.0 m above the surface at 2 meter intervals. The results 
are plotted for each field component using the MATLAB pcolor command. 

The configuration for this example is sketched below. 



02=0-05 mhos/m 


Region 2 



The MATLAB code for this calculation is: 


% Set flags to calculate both E Field and H Field 

EFieldOn = l;HField0n = 1; 

% Frequency = 2 kHz and set dielectric/conductivity for Region 1 & 2 

freqVal = 2000; 
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eDielecl = 10;Conductl = 0.01; 
eDielecZ = 20;Conduct2 = 0.05; 

% Region 1 <--> Air at z = 0 m; Region 2 <--> Region 1 at z = -100 m; 

zLayerl = 0;zLayer2 = -100; 

% Electric dipole moment = 0.5 Amp-m, along y-axis 

IL = 0.5; 

dipMom = IL.*[0,1,0]; 

% Electric dipole moment at x = 0 m, y = 0 m, z = -0.5 m 

% Since the Region 1 <--> Air interface is at z = 0 m (zLayerl = 0 m), 

% this places the dipole o.5 m below the earth 

RSourcePoint = [0,0,-0.5]; 

% Set up a square grid of 30 m x 30 m at a distance of 1.0 m above the 

% Region 1 <--> Air interface at 2 meter intervals on which to 

% evaluate the electric and magnetic fields 

xx = C(-30):2:30).';yy = ((-30):2:30).';zz = 1.0; 

xxTemp = xx*ones(size(yy.'));yyTemp = ones(size(xx))*(yy.'); 

xxEval = xxTemp(:);yyEval = yyTemp(:); 

npts = max(size(xxEval)); 

REvalPoints = zeros(npts,3); 

REvalPoints(:,1) = xxEval;REvalPoints(:,2) = yyEval; 

REvalPoints(:,3) = zz.*ones(npts,1); 

% Do calculation to get E/H Fields at the points on the evaluation grid 

[EField,HField] = getElecDipoleField(EFieldOn,HFieldOn,dipMom, . . . 

RSourcePoint,REvalPoints,freqVal,... 

zLayerl,eDielecl,Conduct1,zLayer2,eDielec2,ConductZ); 


{ MATLAB code to 


do plots } 




Note that it is not necessary to just use points calculated in the x-y plane; any set of 
points is acceptable, as long as all the points in a single call are located in a single layer - 
either Air or Region 1. We selected the set of points to be in the x-y plane for this 
example, just to show one useful configuration. 

The results of this calculation are shown on the following page in the four color 
intensity plots for the three E-Field components and the magnitude of the E-Field along 
with similar plots for the H-Field. We have supplied the complete MATLAB code, 
ExampleHED.m, including the plotting sequences, for this routine on the CD-ROM. 
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Example: Generic Geophex GEM Sensor and a Buried Iron Sphere 

For the second example we examine a calculation for a sensor system similar to the 
GEM sensors developed by Geophex, Ltd. of Raleigh, North Carolina (our 
subcontractor in this effort). The sensor system and measurement configuration are 
sketched below: 


z 



For aM regions: 

b= P 0 

More details on the GEM sensors are given in our Phase 1 Final Report and are, of 
course, available from Geophex. We are using a somewhat generic version of the GEM 
in this calculation and the results will therefore be representative in a general sense, but 
are not meant to be a precise model for any specific version of the GEM sensor. 

In this example a 0.5 m radius iron sphere is buried 1.5 meters below the surface and 
excited by a source coil at one end of the sensor system. This induces a magnetic dipole 
moment on the sphere, which reradiates and is sensed by the receiver coil at the other 
end of the sensor system. The receiver also picks up the direct field from the source coil 
and so the total signal is a sum of the direct field (the field if there sphere were not 
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present) and the perturbed field due to the sphere. Since the GEM is designed to 
subtract out the direct signal, the effective ''measured" signal is the voltage ratio of the 
perturbed signal to the direct signal; this is usually given in parts-per-million (PPM). 

The calculation is done at 1200 Hz; the relative dielectric constant of Region 1 is 10 and 
the conductivity is 0.01 mhos/m; the relative dielectric constant of Region 2 is 50 and 
the conductivity is 0.1 mhos/m. The interface between the Air and Region 1 
corresponds to z = 0 m; the interface between the Region 1 and Region 2 corresponds 
to z = -3 m. We assume the magnetic dipole strength of the source coil is 3 amp-m 2 . The 
sphere is centered at x = 0 m, y = 0 m, and z = -1.5 m in our coordinate system. The z 
component of the perturbed and direct H-Field are evaluated on a square grid 6m x 6m 
a distance of 1.0 m above the surface at ~0.2 meter intervals. The results are plotted for 
each field component using the MATLAB pcolor command. 

The MATLAB routine calls a function getMagDipMomCond. This function is supplied on 
the CD-ROM. The function calculates the induced magnetic dipole moment coefficient 
for a conducting sphere for a given permeability and dielectric constant. The relation 
was derived by Wait ( Geophysics, 16, No. 4, 666-672, October 1951 and Radio Science, Vol. 3 
(New Series), No. 10,1030-1034 October 1968). 

The MATLAB code for this calculation is: 

%***************************************************************************** 
SommerfetdGlobats; 

% Source Loop strength 3 amp*m A 2 and the loop is held to be 

% parallel to the earth, so the dipole moment is alonq the z-axis 

IA_GEM = 3; 

dipMomGEM = IA_GEM.*[0,0,1] ; 

% Source and sensor loops are separated by 2 meters along the y-axis 
LGEM = 2; 

% The initial position of the source is 1 meter (z = 1 meter) above 

% the earth; it is separated by LGEM from sensor along y-axis. 

RSourcePoint0 = [0,-LGEM/2,1]; 

REvalPoints0 = [0,+LGEM/2,1]; 

% Set up measurement grid. 

LL = 3 ;NN = 30; 

nx = NN;xMin = -LL;xMax = +LL;dx = 2*LL/(nx-l); 
xx = xMin:dx:xMax; 
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ny = NN;yMin = -LL;yMax = +LL;dy = 2*LL/(ny-l); 
yy = (yMin:dy:yMax).'; 

freqVal = 1200; 

wVal = 2*pi*freqVal; 

eDielecl = 10;Conductl = 0.01; 
eDielec2 = 50;Conduct2 = 0.1; 

zLayerl = 0;zLayer2 = -3; 

% Set flags to calculate only the H Field 

EFieldOn = 0;FIFieldOn = 1; 

% Get Direct Field at Sensor for each of the measurement position points 
% on the grid 
REvalPoints = zeros(l,3); 

HFieldDirZ = zeros(ny,nx); 

for iy=l:ny 

for ix=l:nx 

gridPoint = [xx(l,ix),yy(iy,1),0]; 

RSourcePoint = RSourcePoint0 + gridPoint; 

REvalPoints = REvalPoints0 + gridPoint; 

[EFieldDirTemp,HFieldDirTemp] = ... 

getMagDipoleField(EFieldOn , FIFieldOn , dipMomGEM, . . . 
RSourcePoint,REvalPoints,freqVal,... 

zLayerl,eDielecl,Conductl,zLayer2,eDielec2,Conduct2); 
HFieldDirZ(iy,ix) = HFieldDirTemp(l,3); 

end 

end 

% Get Field at Underground Sphere for each of the measurement position 

% points on the grid. Sphere located at x=y=0, z=-1.5 m 

REvaUG = [0,0,-1.5]; 

HFieldUGX = zeros(ny,nx); 

HFieldUGY = HFieldUGX;HFieldUGZ = HFieldUGX; 
for iy=l:ny 

for ix=l:nx 

gridPoint = [xx(l,ix),yy(iy,1),0] ; 

RSourcePoint = RSourcePoint0 + gridPoint; 
[EFieldUGTemp,HFieldUGTemp] = ... 

getMagDipoleField(EFieldOn,HFieldOn,dipMomGEM,... 
RSourcePoint,REvaUG,freqVal,... 

zLayerl,eDielecl,Conductl,zLayer2,eDielec2,Conduct2); 
HFieldUGX(iy,ix) = HFieldUGTemp(l,1); 

HFieldUGY(iy,ix) = HFieldUGTempCl,2); 

HFieldUGZ(iy,ix) = HFieldUGTempCl,3); 

end 

end 

% Get Perturbed Field at Sensor for each of the measurement position 

% points on the grid 

epsFreeSpace = 8.85E-12;muFreeSpace = 4E-7*pi; 

epsOut = epsFreeSpace*(eDielecl+i*Conductl/(wVal*epsFreeSpace)); 
muOut = muFreeSpace; 

siglron = 1.03E7;muIron = 3*muFreeSpace; 

epsSph = epsFreeSpace*(0.0+i*sigIron/(wVal*epsFreeSpace)); 
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muSph = mulron; 

ax = 0.5;ay = 0.5;az = 0.5; 

[inducedDipMomCoef] = ... 

getMagDipMomCond(freqVal,epsOut,muOut,epsSph,muSph,ax,ay,az); 

RSourceUG = REvaUG; 

HFieldPertbZ = zeros(ny,nx); 
for iy=l:ny 

for ix=l:nx 

gridPoint = [xx(l,ix),yy(iy,1),0]; 

REvalPoints = REvalPoints0 + gridPoint; 

HFieldUG = [HFieldUGX(iy ,ix),HFieldUGY(iy,ix),HFieldUGZ(iy,ix)]; 
inducedDipMom = inducedDipMomCoef. *FIFieldUG; 
[EFieldPertbTemp,FIFieldPertbTemp] = ... 

getMagDipoleField(EFieldOn,HFieldOn,inducedDipMom, . . . 
RSourceUG,REvalPoints,freqVai,... 

zLayerl,eDielect,Conduct1,zLayer2,eDielec2,Conduct2); 
HFieldPertbZ(iy,ix) = FIFieldPertbTemp(l, 3) ; 

end 

end 


% Find absolute Perturbed/Direct (in PPM) for H Field for each of the 
% measurement position points on the grid 
HFieldRatioPPM = 1.0E6.*abs(HFieldPertbZ./HFieldDirZ) ; 

% Plot Perturbed H Field in PPM for each of the measurement position 

% points on the grid 

figure 

pcolor(xx,yy,HFieldRatioPPM) 
shading interp 
xlabel('x - meters') 
ylabelC'y - meters') 
colorbar 

titleStrLl = sprintf(' Hz(perturb)/Hz(direct) [PPM]\n'); 
titleStrL2 = 'A 0.5 m Radius Iron Sphere at 1.5 m below surface'; 
titleStr = [titleStrLl,titleStrL2]; 
title(titleStr) ; 
axis image 

%***************************************************************************** 

We have supplied this MATLAB code, LoopSensorExample.m, along with 
getMagDipMomCond.m on the CD-ROM. 


The results of this calculation are shown below as a color intensity plot of the ratio 

| Perturb j |_J Direct | 

for the grid in PPM; the results have been smoothed using the MATLAB shading 
interp command. 
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Menu Interface Routines 

This section describes the MATLAB interface that allows a user to define the sources, 
receivers and buried structures. The user can then setup appropriate linkages to any 
computations to configure an easy to use menu driven interface. The menu interface 
software has been developed so that linkages to any user defined MATLAB code 
(including the Sommerfeld routines or the Hill/King buried wire routines) may be 
easily implemented by an experienced MATLAB user. 

Objects 

1. There are objects called "sources" and "receivers" that have locations that vary; i.e., 
they can move. These locations are called "measurement points" or "time steps". 
However, the order of the measurement points is not important; they are not time 
related. There may be any number, up to some maximum, of measurement points. 
There must be at least one measurement point. 

There may be any number, up to some maximum, of the objects. For example, you 
may have 1 source and 3 receivers. All sources and receivers have the same number of 
measurement points. There must be at least one source and one receiver. 

2. There are also objects known as "buried structures" that have locations that do not 
vary; i.e., they cannot move. Again there may be any number, up to some maximum, 
of these objects. And, there may be no (0) buried structures. 

3. A location is specified in Cartesian coordinates as a point, (x, y, z), in meters 
referenced from an arbitrary origin, (0, 0, 0). A measurement point is also specified in 
meters, (x, y, z), from the same arbitrary origin. 

4. All sources have a strength that is given in units appropriate to the type of the source 
There are four types of sources: (1) loop, (2) electric dipole, (3) magnetic dipole and (4) 
straight wire. 

a. A loop source is defined by its radius, in meters, and its orientation, the plane that is 
perpendicular to a specified unit vector located at its center point. The unit of strength 
of a loop source is "Amp". The location of the center of the loop will define the location 
of the loop. The loop intended here has finite size with respect to other dimensions of 
the problem, e.g. layer depths, distances to objects, etc. 
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b. An electric dipole source is located at a single point and is defined by its orientation, a 
specified unit vector at that point. The unit of strength of an electric dipole source is 
"Amp-m". It is an infinitesimal-length wire fed at its center, in terms of its radiated fields. 

c. A magnetic dipole source is also located at a single point and defined by its 
orientation, a specified unit vector at that point. The unit of strength of a magnetic 
dipole source is "Amp-m 2 ". It is an infinitesimally small loop, whose axis coincides with 
the perpendicular to the loop center, with respect to its radiated fields. 

d. A straight wire source is defined by its length and direction, both given by a 
specified vector with center at the midpoint of the wire. The location of the center of 
this vector will define the location of the wire. The unit of strength of a wire source is 
"Volt". If the wire length is short (approaching infinitesimal), the wire is usually better 
or more conveniently modeled by the electric dipole. The feed of the wire is assumed to 
be at the center (this may be changed by a MATLAB modification to the menu interface 
program). 

5. There are four types of receivers, the same as the sources. However, receivers do not 
have an associated strength. All other properties are the same as the four 
corresponding types of sources. 

6. There are two types of buried structures: (1) ellipsoid and (2) wire. 

a. An ellipsoid buried structure is defined by its 3 semi-axes, ax, ay and az; its 
orientation, which may be specified by two unit vectors (one corresponding to its local 
z-axis and one corresponding to its local x-axis) or three rotations (angular rotations 
about each of its local axes); and its electrical properties which are its relative 
permeability, relative dielectric constant or permittivity, and whether or not it is a 
perfect conductor. The center of the ellipsoid will define its location. 

b. A straight wire buried structure is defined by its length and direction, both given by 
a specified vector with center at the midpoint of the wire. The center of this vector will 
define its location. 

Global Variables 

The above information is communicated to the rest of the program via the following 
global variables, found in the m-file "SRLGlobals.m". 
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Sources / Receivers 

nMeasure - the number of measurement points. 
nSource - the number of sources. 
nReceiver - the number of receivers. 

xSource, ySource, zSource - arrays (nMeasure,nSource). 

These arrays contain the measurement locations in meters of all sources. 

sourceType - array (l,nSource). 

This array contains the type of each source, 1 == loop, 2 == electric dipole, 
3 == magnetic dipole and 4 == wire. 

sourceStrength - array (l,nSource). 

This array contains the strength of each source. 

sourceProperties - array (nSource,4). 

This array contains the defining properties of each source. There is one 
row per source and the contents of each row varies with the type of 
source as follows: 

loop - columns 1-3 == the orientation unit vector, 

column 4 == radius of the loop in meters; 
dipoles - columns 1-3 == the orientation unit vector, 
column 4 == not used; 

wire - columns 1-3 == the orientation and length vector, 

column 4 == not used; 

xReceiver, yReceiver, zReceiver - arrays(nMeasure,nReceiver). 

These arrays contain the measurement locations in meters of all receivers. 

receiverType - array (l,nReceiver). 

This array is similar to the sourceType array. 

receiverProperties - array (nReceiver,4). 

This array is similar to the sourceProperties array. 
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Additional source / receiver variables 

The following variables are created if the user decides to have measurement locations 
computed from a grid. 


xGrid, yGrid, zGrid - arrays (4). 

These arrays contain the definition of the computed grid. The four 
elements of each array define the minimum value in meters, the 
maximum value in meters, the number of grid points and the distance 
between grid points in meters - for each axis. 

xOSource, yOSource, zOSource - arrays(nSource). 

These arrays contain the location of each source relative to the origin of 
the grid. 

xOReceiver, yOReceiver, zOReceiver - arrays (nReceiver). 

These arrays contain the location of each receiver relative to the origin of 
the grid. 

Buried Structures (Ellipsoids) 

nElip - the number of ellipsoids defined. 

xElip, yElip, zElip - arrays (l,nElip). 

These arrays give the locations, in meters, of the centers of the various 
buried ellipsoids. 

ixElip, iyElip, izElip - arrays (l,nElip). 
jxElip, jyElip, jzElip - arrays (l,nElip). 
kxElip, kyElip, kzElip - arrays (l,nElip). 

These arrays describe three unit vectors for each ellipsoid. The unit 
vectors give the orientation of each ellipsoid. The unit vector (jxElip, 
jyElip, jzElip) is computed by the program given the other two. 

dxElip, dyElip, dzElip - arrays (l,nElip). 

These arrays give the lengths, in meters, of each of the three semi-axes of 
each ellipsoid. 
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pcElip - array (l,nElip). 

This array provides a "perfect conductor" flag for each ellipsoid. (== 0 
means the ellipsoid is not a perfect conductor). 

epsElip - complex array (l,nElip). 

This array provides the relative complex dielectric constant for each 
ellipsoid. 

muElip - complex array (l,nElip). 

This array provides the relative complex permeability of each ellipsoid. 

Buried Structures (Wires) 

nWire - the number of wires defined. 

xWire, yWire, zWire - arrays (l,nWire). 

These arrays give the locations, in meters, of the centers of the various 
buried wires. 

vxWire, vyWire, vzWire - arrays (l,nWire). 

These arrays describe a vector for each wire. The center of the vector is 
located at the center of its corresponding wire. And the direction of the 
vector is the direction of the wire. And the length of the vector is the 
length of the wire (meters). 

Private Global Variables 

The source / receiver and buried structures information is organized into arrays of 
structures for the use of the input portion of the program. These arrays of structures 
are defined within a "private" global area. (See the m-file "SRLPrivates.m".) Here are 
descriptions of these private global variables. 

Sources / Receivers 

SRLmethod - the method of input, 0 == by computation, 1 == by file. 

SRLpathname - the path and name of the input file. 
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This variable is used if the input is by file. 

SRLgridX, SRLgridY, SRLgridZ 

These arrays hold the same information as the variables xGrid, yGrid and 
zGrid described above. 

SRLSn - the number of sources. 

SRLSmax - the maximum allowed number of sources. 

SRLSrc - An array of structures (a cell array). 

Each structure within this array contains the relative measurement 
location and defining properties of one source. There is one 
structure per source. The cell array is preallocated with the 
MATLAB statement "SRLSrc = cell(SRLSmax,l);". 

Here are the elements of each structure: 

.type = integer - the type of source, 

1 == loop 

2 == electric dipole 

3 == magnetic dipole 

4 == wire 

.location = array (1,3) 

The location (x,y,z - meters) of this source relative to 
the current measurement point 
.vector = array (1,3) the orientation (and/or length) vector 
.radius = the radius of a loop source. 

.strength = complex - the source strength; 

Specific quantities within the array of structures are referred to 
with notation in the form of "SRLSrc{5}.location(2)", which refers to 
the y-coordinate of the relative location of the 5th source. 

SRLRn - the number of receivers. 

SRLRmax - the maximum allowed number of receivers. 
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SRLRcv - this is a cell array of structures that define the locations and 
of the receivers. It is defined similar to the above cell array, 
difference is that the value of the "strength" member of the 
within the array will be set to zero. 

Buried Structures 

BSDpathname - the path and name of the input file. 

This variable is used if the input is by file. 

BSDEn - the number of ellipsoids. 

BSDEmax - the maximum allowed number of ellipsoids. 

BSDEllip - cell array of ellipsoid buried structures 

One structure per ellipsoid defined as follows: 

.location = array (1,3)/ the location of the center of the ellipsoid in 
Cartesian coordinates (x,y,z - meters) 

.axes = array (1,3), the lengths of the semi-axis, x, y and z, in 
meters 

,xunit= array (1,3), unit vector describing the orientation of the 
ellipsoid's x-axis relative to its local origin 
,yunit= array (1,3), unit vector describing the orientation of the 
ellipsoid's y-axis (this vector is calculated) 

,zunit= array (1,3), unit vector describing the orientation of the 
ellipsoid's z-axis relative to its local origin 
.angles = array (1,3), alternate orientation expressed as angles of 
rotation about the respective axes (degrees), (the x,y,z-unit 
vectors are computed from this), 

.orient = integer, orientation flag, l==unit vectors, 0==angles, 

.pc = integer, perfect conductor flag, l==yes, 0==no, 

.eps= complex, relative dielectric constant, 

.mu = complex, relative permeability constant. 

BSDWn - the number of wires. 


properties 
The only 
structures 
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BSDWmax - the maximum allowed number of wires. 


BSDWire - cell array of wire buried structures 

One structure per wire defined as follows: 

.location = array(l / 3) / the location of the midpoint of the wire in 
Cartesian coordinates (x,y,z - meters) 

.vector = array(l / 3), a vector with its center at the midpoint of the 
wire. This vector describes length (meters) and direction. 


Operation 

To input source and receiver locations invoke the MATLAB function "SRLocations". To 
input the buried structures definitions invoke the MATLAB function "BSDefinition". 

Each of these functions will present a dialog box that allows the user to select between 
input by computation or input by file. 

Input via File 

If you intend to select the input by file method you must have on the disk a file that is 
readable via the MATLAB "load filename" statement. 

If your file input is for the sources and receivers then your input file must have all of the 
source and receivers variables described in the sections "Sources/Receivers" and 
"Additional source/receiver variables" in the "Globals" section above. 

If your file input is for the buried structures then your input file must have all of the 
buried structures' variables described in the sections "Buried Structures (Ellipsoids)" and 
"Buried Structures (Wires)" in the "Globals" section above. 

Source / Receiver Input via Computation 

When you input your source/receiver definitions via computation you must define the 
grid of measurement locations, define at least one source and define at least one 
receiver. At the conclusion of your input the definitions will be converted and copied to 
the global variables described in the "Globals" section above. These values will also be 
written to the file "SRLSave.mat" in a form that is readable via the MATLAB "load" 
statement. 
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Each time the "SRLocations" function is invoked it will look for the "SRLSave.mat" file 
and read it and present those values to you when you select input by computation. If 
you want to start from scratch discard the file "SRLSave.mat". 

The dialog boxes presented during the input by computation are self explanatory. Each 
box represents one source or one receiver. You navigate from source to source or 
receiver to receiver via the buttons within the box. You can choose to add or delete 
items. When you delete an item it is moved to the end of the list and the list is 
shortened by one. Thus, when you subsequently add a new item, the one you deleted 
appears at the end of the list for user information. 

Buried Structures Input via Computation 

When you input your buried structures definitions via computation the process is 
similar to that described above. At the conclusion of your input the definitions will be 
converted and copied to the global variables described in the "Globals" section above. 
And these values will also be written to the file "BSDSave.mat" in a form that is readable 
via the MATLAB "load" statement. 

Each time the "BSDefinition" function is invoked it will look for the "BSDSave.mat" file 
and read it and present those values to you when you select input by computation. 

And, if you want to start from scratch then discard the file "BSDSave.mat". 

The dialog boxes presented for buried structure definitions work the same as those for 
source / receiver definitions. There is one box for each item. You may have zero items, 
in which case all entries in the box are disabled except the "add" button. 
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Hill/King Buried Wire Routines 

For a number of applications it is useful to model long insulated conductor in a lossy 
medium with arbitrary termination on either end. We have written MATLAB routines 
that use David A. Hill's paper (expanding on previous work of R.W.P. King - "Magnetic 
Dipole Excitation of an Insulated Conductor of Finite Length", IEEE Trans. Geosci. 
Remote Sensing, vol. 28, pp. 289-294,1990) to do these calculations. The Hill/King result 
allows current on an insulated conductor in a lossy medium with arbitrary termination 
on either end to be obtained from a defined, imposed, external electric field. 

The Hill/King model assumes that the wire is buried in a single infinite medium; we 
are interested in multilayered applications - so there is a issue of how the current is 
affected by getting near a medium interface; we have used NEC-3 to address this issue 
(discussed in detail below). It is our opinion that the single infinite medium 
approximation for the calculation of the Hill / King current will work well for multilayer 
applications. Included on the CD-ROM are examples using this code; the examples are 
taken from Hill's paper reproducing Figures 6-10 in the paper. 

To determine the effect of a medium interface on a current distribution we used NEC-3 
to look at the current distribution as a voltage excited horizontal buried wire moves 
closer to a medium interface. We used 3 kHz and average earth (s=10, o=0.01 mhos/m). 
We looked at two lengths of wire: one a little under 2.5 k medium in length (1000 m) and 
one a little under 0.5 A, medium in length (200 m). The wire was fed at the center with 1 Volt. 
The depths used were 2000 m, 200 m, 20 m, 2 m, and 0.2 m. The results of these 
calculations show that there is a small change in the magnitude of the current and a 
slight increase in the phase as the wire approaches the interface. 
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Current Magnitude for Center Fed, Horizontal, Buried Wire at Various Depths Using NEC-3 at 3 kHz 
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Current Phase for Center Fed, Horizontal, Buried Wire at Various Depths Using NEC-3 at 3 kHz 
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Current Magnitude for Center Fed, Horizontal, Buried Wire at Various Depths Using NEC-3 at 3 kHz 
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Current Phase for Center Fed, Horizontal, Buried Wire at Various Depths Using NEC-3 at 3 kHz 
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VED in Region 1 

Located at z = 0 - Interface 1-2 at z = -d l - Interface 2-3 at z = -d 2 


Fields in Region 1: (z > -d,) 
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VMD in Region 1 

Located at z = 0 - Interface 1-2 at z = -d l - Interface 2-3 at z = -d 2 

Fields in Region 1: (z > -d,) 
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VED in Region 2 

Located at z = -Zq (-d 2 < -z 0 < -d,) - Interface 1-2 at z = -d l 


Fields in Region 1: (z > -d,) 
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Fields in Region 2 above the source : (-d, > z > -z 0 ) 
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Interface 2-3 at z = -d 2 
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VMD in Region 2 

Located at z = -4 (-d 2 < -z 0 < -dj) - Interface 1-2 at z = -d. 


Fields in Region 1: (z > -d,) 
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Fields in Region 2 above the source : (-d, > z > -z 0 ) 
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HED in Region 1 

Located at z = 0 

Fields in Region 1: (z > -d x ) 


Interface 1-2 at z = -d, 


Interface 2-3 at z = -d. 
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Fields in Region 2 : (-dj > z > -d 2 ) 
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HMD in Region 1 

Located at z = 0 

Fields in Region 1: (z > -d,) 
IA 


Interface 1-2 at z = -d, 


Interface 2-3 at z = -d. 
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Fields in Region 2: (-dj > z > -d 2 ) 


H 
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sin^T^/^ p —^+ fl“A ( ” 2 ‘ W ) 

4jt £ ^ 


j^HMD[ 1] _ . Pi^ 4 


2 p 


/-COS 


4mi 2 


oo 

<t>fdk p A n: k 2 M f , J Q (k p p)~ 


J '{ k pP)\ 


k p P j 


_e~ ik 2z z + ftTE e ik 2z (z+2<k )| 


~^ 2 /' |M —COS(pcdk A™ -i-pY 

%2 ^2 { K V 


e ~ik 2z z + ftTM e ik 2z (z+2d 2 ) 


H™ m = i ^ sin</>J dk 


20 


4jrpu 2 


^J{k p p\e^-^^) 


+ik■ 


; / 2 ,U|M £ | 

4jrp 2 £ 2 


sin 


4“t f 


0 \ 


2 ,(*,p) 


V / 


p-ikizZ _j_ ft™ g ik i= (z +2d i) 


C 01 ' 1 - + R!je 

4JTO *o 

00 

l c/jj' dk 


ik 2z (z+2d 2 ) 


4jip 

. CO Id J A £, 


sin 4. 

4jr ^2 J 0 








+ R™<? 2 ‘ 


ik 2 _(z+2d 2 ) 


4jt 
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00 
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ELF/VLF SBIR Phase II Final Report 


CODAR Ocean Sensors, Ftd. Page 34 of 42 



f,TM jjTM 

‘ 2\ £*2 


HED in Region 2 

Located at z = -4 - Interface 1-2 at z = - 

Fields in Region 1: (z > -d x ) 

= i _U_ COS( pr dk k 2 j( k p \ e 

4jtw£, 4 ^ p v ^ ' 

11 “ t 2 

Ait u, J k 2 , 

——cos^f dk -!- J^k py^^'V^'f™B™ 

4jt0)£ i p P , *» A 

J1 00 
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Interface 2-3 at z = 
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Fields in 


pH£D[ 2 ] • II 

C 2 z 


source: (-dj > z> -z 0 ) 


—^_ili cos 0 r(ik -^l y (t p], 
4 tt P Pl J p k 2z p > 

Region 2 above the 
II 00 

T1 °° k 2 
i—sin^fdkp -^J t (k p p)B: 
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r^HEDf 21 II 

, L J = - 
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Fields in Region 2 below the source: (-4 > z > -d 2 ) 
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K^ 2] = _ ~7~~ sin <pCdk f) J l [k p p ^C™{e~ ik2zZ + e^^R™) 

4jrp J Q 
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HMD in Region 2 

Located at z = -Zg - Interface 1-2 at z = -d, - Interface 2-3 at z = -d 2 


Fields in Region 1: (z > -d,) 
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Fields in Region 2 above the source: (-dj > z > -z 0 ) 
K MD[2] = ~~~^ n( Pjd k p -^-J^kpPjM™^ 
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+ sin (pjdk p M™ k p 

o 


4jt 


■4>M 




V 


£ p p 


git 2 Z Z _ ^-ik 2z (z+2d l )j^TM 


ELF/VLF SBIR Phase II Final Report 


CODAR Ocean Sensors, Ltd. Page 38 of 42 



E T [2] = 1 ° )! ^ A cos <pfdk p M T i E k p J 0 (k p p)~ 

+ 7^1 c 0 s<Pf dkp M™J t (k pP je lk ^ - e ~ ik ^ +2d ')R™ 
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Fields in Region 2 below the source: (-4 > z > -d 2 ) 
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= -—sin ^dk p k 2i J x [k p p)Nl E [-e~ ,k ^ z + e 


jk lz (z+2d 2 ) nTE 
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Appendix B - Contents of CD-ROM 

On the CD-ROM, named ELF_VLF delivered with this report are three 
folders/ subdirectories: 

Hill_King, Menujnterface, Sommerfeld_Routines 

In addition to these folders/subdirectories is a PDF file version of this Report 

ELF_Report.PDF” 

Below are listed the files contained in each of the folders/subdirectories 

Hill_King: 

HillFiglO.m 

HillFig6.m 

HillFig7.m 

HillFig8.m 

HillFig9.m 

KingsGreenFtn.m 

Menulnterface: 

BSDefinition.m 

BSDefinitionBox.m 

BSDefinitionBox.mat 

BSDEIlipsoid.m 

BSDEIlipsoidBox.m 

BSDEIlipsoidBox.mat 

BSDWireStruct.m 

BSDWireStructBox.m 

BSDWireStructBox.mat 

ellipplot.m 

SRLGIobals.m 

SRLGrid.m 

SRLGridBox.m 

SRLGridBox.mat 

SRLObject.m 

SRLObjectBox.m 

SRLObjectBox.mat 

SRLocations.m 

SRLocationsBox.m 

SRLocationsBox.mat 

SRLPrivates.m 

Sommerfeld_Routines: 

EHEDRegion1_1 .m 
EFIEDRegion1_2.m 
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EHEDRegion2_1 .m 

EHEDRegion2_2.m 

EHMDRegion1_1 .m 

EHMDRegion1_2.m 

EHMDRegion2_1 .m 

EHMDRegion2_2.m 

EVEDRegion1_1 .m 

EVEDRegion1_2.m 

EVEDRegion2_1 .m 

EVEDRegion2_2.m 

EVMDRegion1_1 .m 

EVMDRegion1_2.m 

EVMDRegion2_1 .m 

EVMDRegion2_2.m 

ExampleHED.m 

fasti nt.m 

fastlntSklz.m 

getElecDipoleField.m 

getkrlnterp.m 

getMagDipMomCond.m 

getMagDipoleField.m 

FIHEDRegion1_1 .m 

HHEDRegion1_2.m 

HHEDRegion2_1 .m 

HHEDRegion2_2.m 

HHMDRegion1_1 .m 

HHMDRegion1_2.m 

HHMDRegion2_1 .m 

HHMDRegion2_2.m 

FIVEDRegion1_1 .m 

HVEDRegion1_2.m 

HVEDRegion2_1 .m 

HVEDRegion2_2.m 

HVMDRegion1_1 .m 

HVMDRegion1_2.m 

HVMDRegion2_1 .m 

HVMDRegion2_2.m 

interpl COS.m 

LoopSensorExample.m 

psiN.m 

SommerfeldGlobals.m 

Sommerfeldlnit.m 
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